% May 23, 2025
% Duffy, Lebeau, Puzzello

% This file replicates Table 10 and Figure C1.

%% Set directory HERE
%Example: cd 'C:\Users\Lenovo\Desktop\replication_DLP\utility'
cd 'C:\...\replication_DLP\utility'

%% Load data
opts = spreadsheetImportOptions("NumVariables", 3);
opts.Sheet = "Sheet1";
opts.DataRange = "A2:C102";
opts.VariableNames = ["Year", "InterestRate", "MoneyDemand"];
opts.VariableTypes = ["double", "double", "double"];
dataCraigRocheteau = readtable("data_CraigRocheteau.xlsx", opts, "UseExcel", false);
clear opts

data=table2array(dataCraigRocheteau);
interest=data(:,2)/100;
money=data(:,3);

%% Estimate benchmark: model with theta=0.5 and Kalai
theta=0.5;
kalai=@(x,y) modelkalai(y,x(2),theta,x(1)); 
x0 = [0.01,0.01];
lb = [0,0.1];
ub = [500,1];
solkalaibench = lsqcurvefit(kalai,x0,interest,money, lb, ub);

%% Figure C1
interestgrid=0.0001:0.001:0.15;
x=0;
for r=interestgrid
    x=x+1;
    xvalues(x)=modelkalai(r,solkalaibench(2),theta, solkalaibench(1));
end 
a=modelkalai(0.03,solkalaibench(2),theta, solkalaibench(1));
b=modelkalai(0.13,solkalaibench(2),theta, solkalaibench(1));
set(0,'DefaultFigureColor',[1 1 1]);

figure
plot(money,interest,'o','LineWidth',1.1,'Color',[0,0,0])
hold on
plot(xvalues,interestgrid,'LineWidth',2)
hold on
plot(0.05:0.0001:a,0.03*ones(length(0.05:0.0001:a),1),'--','LineWidth',2,'Color',[0,0,0]+0.3)
hold on
plot(0.05:0.0001:b,0.13*ones(length(0.05:0.0001:b),1),'--','LineWidth',2,'Color',[0,0,0]+0.3)
xlabel("Money demand")
ylabel("Nominal interest rate")
yticks([0.01:0.02:0.15])
legend('Data','Fitted Curve')
set(gca,'Color','white','FontName','Garamond','LineWidth',1.2,'FontWeight','bold','FontSize',12);

%% Table 10

% Compute benchmark outcomes
[deltabench, qzerobench, qbasebench, qrbench] = welfarecostkalai(0.13,solkalaibench(2),0.5,solkalaibench(1)); 

% Compute cost for different theta with Kalai,reestimating the As and sigmas
solkalai33 = lsqcurvefit(@(x,y) modelkalai(y,x(2),0.33,x(1)),x0,interest,money,lb,ub);
solkalai66 = lsqcurvefit(@(x,y) modelkalai(y,x(2),0.66,x(1)),x0,interest,money,lb,ub);
solkalai100 = lsqcurvefit(@(x,y) modelkalai(y,x(2),1,x(1)),x0,interest,money,lb,ub);
A = [solkalai33(1) solkalaibench(1) solkalai66(1) solkalai100(1)];
sigma = [solkalai33(2) solkalaibench(2) solkalai66(2) solkalai100(2)];

[deltakalai33, qzerokalai33, qbasekalai33, qrkalai33] = welfarecostkalai(0.13,solkalai33(2),0.33,solkalai33(1)); 
[deltakalai66, qzerokalai66, qbasekalai66, qrkalai66] = welfarecostkalai(0.13,solkalai66(2),0.66,solkalai66(1)); 
[deltakalai100, qzerokalai100, qbasekalai100, qrkalai100] = welfarecostkalai(0.13,solkalai100(2),1,solkalai100(1));  
result=100*[deltakalai33 deltabench deltakalai66 deltakalai100]
qzero=[qzerokalai33 qzerobench qzerokalai66 qzerokalai100] % q traded at Friedman rule (optimal)
qbase=[qbasekalai33 qbasebench qbasekalai66 qbasekalai100] % q traded at r=3%
qr=[qrkalai33 qrbench qrkalai66 qrkalai100] % q traded at r=13%

% Compute cost of inflation for different theta with Nash and A fitted for theta=0.5 (Nash)
theta=0.5;
nash=@(x,y) modelnash(y,x(2),theta,x(1)); 
x0 = [0.01,0.01];
lb = [0,0.1];
ub = [500,1];
solnashbench = lsqcurvefit(nash,x0,interest,money);

[deltabenchn, qzerobenchn, qbasebenchn, qrbenchn] = welfarecostnash(0.13,solnashbench(2),0.5,solnashbench(1)); 

% Compute cost for different theta with Nash, reestimating the As and sigmas
solnash33 = lsqcurvefit(@(x,y) modelnash(y,x(2),0.33,x(1)),x0,interest,money);
solnash66 = lsqcurvefit(@(x,y) modelnash(y,x(2),0.66,x(1)),x0,interest,money);
solnash100 = lsqcurvefit(@(x,y) modelnash(y,x(2),1,x(1)),x0,interest,money);
A2 = [solnash33(1) solnashbench(1) solnash66(1) solnash100(1)]; 
sigma2 = [solnash33(2) solnashbench(2) solnash66(2) solnash100(2)]; 

[deltanash33, qzeronash33, qbasenash33, qrnash33] = welfarecostnash(0.13,solnash33(2),0.33,solnash33(1)); 
[deltanash66, qzeronash66, qbasenash66, qrnash66] = welfarecostnash(0.13,solnash66(2),0.66,solnash66(1)); 
[deltanash100, qzeronash100, qbasenash100, qrnash100] = welfarecostnash(0.13,solnash100(2),1,solnash100(1));  
result2=100*[deltanash33 deltabenchn deltanash66 deltanash100]
qzero2=[qzeronash33 qzerobenchn qzeronash66 qzeronash100] % q traded at Friedman rule (optimal)
qbase2=[qbasenash33 qbasebenchn qbasenash66 qbasenash100] % q traded at r=3%
qr2=[qrnash33 qrbenchn qrnash66 qrnash100] % q traded at r=13%